Tutorial
LiDAR Tutorial using R
Introduction
Packages
library(lidR)# for working with LAZ files
library(sf) # for working spatial class
library(raster)# for working with raster
library(rayshader) # for 3D viz
library(rgl) # for interactive plotsArea of Interest
motte_point = st_read("data/canmore_point.shp", quiet = TRUE)
motte_buffer = st_buffer(motte_point,dist = 50)LAZ
Next, I am reading the LAZ files and clipping the point cloud to the extent of immidiate area around the motte.
#read laz files
las = readLAS("data/NX6055_4PPM_LAS_PHASE3.laz")
motte = clip_roi(las, motte_buffer)
#plot lidar point cloud
plot(motte, bg = "white")
rglwidget()grab and rotate the model !
Digital Elevation Model
In this step I am running standard algorithms from LiDR package to compute Digital Surface Model (DSM) and Digitial Terrain Model (DTM).
rgl.clear()
# create dsm and dtm
dsm = grid_canopy(motte, 0.5, pitfree())
# assign coordinate system
crs(dsm) = CRS("+init=epsg:27700")
# create dtm
dtm = grid_terrain(motte, 0.5, algorithm = knnidw(k = 6L, p = 2))
# addign coordinate system
crs(dtm) = CRS("+init=epsg:27700")
par(mfrow = c(1,2))
plot(dtm, main = "DTM", col = height.colors(50))
plot(dsm, main = "DSM",col = height.colors(50))
Hillshade
# dsm
slope_dsm = terrain(dsm, opt = 'slope')
aspect_dsm = terrain(dsm, opt = 'aspect')
hill_dsm = hillShade(slope_dsm, aspect_dsm, angle = 40, direction = 270)
# dtm
slope_dtm = terrain(dtm, opt = 'slope')
aspect_dtm = terrain(dtm, opt = 'aspect')
hill_dtm = hillShade(slope_dtm, aspect_dtm, angle = 5, direction = 270)
#plot
par(mfrow = c(1,2))
plot(hill_dtm, main = "DTM Hilshade", col = grey.colors(100, start = 0, end = 1),
legend = FALSE)
plot(hill_dsm, main = "DSM Hillshade", col = grey.colors(100, start = 0, end = 1))
3D Model
#And convert it to a matrix:
elmat = raster_to_matrix(dtm)
elmat %>%
sphere_shade(texture = "imhof1") %>%
add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
max_darken = 1) %>%
add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
phi = 40, theta = 180, zoom = 0.8, fov = 1)
rglwidget()grab and rotate the model !
Another example …
And here is the code
library(rayshader)
library(raster)
library(magick)
library(gifski)
dem = raster("data/tantallon.tif")
n_frames <- 41
elmat = raster_to_matrix(dem)
zscale <- 0.5
waterdepthvalues <- c(0:20, seq(19,0,-1)) / 2
phi_values = seq(20,60)
waterdepthvalues * zscale
length(waterdepthvalues)
img_frames <- paste0("drain", seq_len(n_frames), ".png")
for (i in seq_len(n_frames)) {
elmat %>%
sphere_shade(texture = "imhof1") %>%
add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
max_darken = 1) %>%
add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
water = TRUE, watercolor = "imhof3", wateralpha = 0.5,
waterlinecolor = "#ffffff", waterlinealpha = 0.5,
waterdepth = waterdepthvalues[i],
phi = 50, theta = 45, zoom = 0.8, fov = 1)
render_snapshot(filename = img_frames[i],
title_text = "Tantallon Castle",
title_font = "Helvetica",
vignette = TRUE,
title_size = 50,
title_color = "black")
rgl::clear3d()
}
magick::image_write_gif(magick::image_read(img_frames),
path = "tantallon.gif",
delay = 6/n_frames)sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rgl_0.100.54 rayshader_0.15.1 sf_1.0-2 lidR_3.1.4
## [5] raster_3.1-5 sp_1.4-5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.7 lattice_0.20-41 png_0.1-7
## [4] prettyunits_1.1.1 class_7.3-17 digest_0.6.27
## [7] foreach_1.5.0 utf8_1.2.2 mime_0.11
## [10] R6_2.5.1 evaluate_0.14 e1071_1.7-4
## [13] highr_0.9 pillar_1.6.2 rlang_0.4.11
## [16] progress_1.2.2 lazyeval_0.2.2 data.table_1.14.0
## [19] miniUI_0.1.1.1 jquerylib_0.1.4 magick_2.3
## [22] rmarkdown_2.10.1 rgdal_1.5-18 webshot_0.5.2
## [25] stringr_1.4.0.9000 htmlwidgets_1.5.1 munsell_0.5.0
## [28] shiny_1.4.0.2 compiler_4.0.2 httpuv_1.5.4
## [31] xfun_0.25 pkgconfig_2.0.3 rgeos_0.5-2
## [34] htmltools_0.5.1.1 tidyselect_1.1.0 tibble_3.1.4
## [37] codetools_0.2-16 fansi_0.5.0 crayon_1.4.1
## [40] dplyr_1.0.2 later_1.1.0.1 grid_4.0.2
## [43] jsonlite_1.7.2 xtable_1.8-4 lifecycle_1.0.0
## [46] DBI_1.1.0 magrittr_2.0.1 scales_1.1.1
## [49] units_0.6-7 KernSmooth_2.23-17 stringi_1.7.3
## [52] promises_1.1.1 doParallel_1.0.15 bslib_0.2.5.1
## [55] ellipsis_0.3.2 generics_0.1.0 vctrs_0.3.8
## [58] iterators_1.0.12 tools_4.0.2 manipulateWidget_0.10.1
## [61] glue_1.4.2 purrr_0.3.4 hms_0.5.3
## [64] crosstalk_1.1.0.1 parallel_4.0.2 fastmap_1.0.1
## [67] yaml_2.2.1 colorspace_2.0-2 rlas_1.3.5
## [70] classInt_0.4-3 knitr_1.33 sass_0.4.0